R programming: generics

Published

March 21, 2025

M1 MIDS/MFA/LOGOS

Université Paris Cité

Année 2024

Course Homepage

Moodle

Code
stopifnot(
  require(Hmisc),
  require(skimr),
  require(patchwork),
  require(ggforce),
  require(glue),
  require(ggfortify),
  require(broom),
  require(tidyverse)
)

tidymodels::tidymodels_prefer(quiet = TRUE)

old_theme <-theme_set(
  theme_minimal(base_size=9, 
                base_family = "Helvetica")
)
Code
gc <- options(ggplot2.discrete.colour="viridis")
gc <- options(ggplot2.discrete.fill="viridis")
gc <- options(ggplot2.continuous.fill="viridis")
gc <- options(ggplot2.continuous.colour="viridis")

Objectives

Code
stopifnot(
  require(rlang),
  require(lobstr),
  require(sloop),
  require(devtools),
  require(usethis),
  require(testthat),
  require(generics)
)
Loading required package: rlang
Loading required package: lobstr
Loading required package: sloop
Loading required package: devtools
Loading required package: usethis
Loading required package: testthat
Loading required package: generics

Generics and S3 classes

Let us first create an instance of class lm.

Code
lm0  <- lm(Gas ~ Insul * Temp, MASS::whiteside)
Question
  • What does function class() do?
  • Is it possible to belong to type list and to class lm simultaneously?
  • In R what is an attribute?
  • How do we set and get attributes?
  • What does function inherits() do?
Solution
Code
class(lm0)
[1] "lm"
Code
is.list(lm0)
[1] TRUE
Code
attributes(lm0)
$names
 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "contrasts"     "xlevels"       "call"          "terms"        
[13] "model"        

$class
[1] "lm"
Code
inherits(lm0,"list")
[1] FALSE
Code
inherits(lm0, "lm")
[1] TRUE
Question

Load package sloop.

  • What does sloop::otype() do? Apply it to an object of class lm.
  • What happens when we first unclass() the object?
Solution
Code
sloop::otype(lm0)
[1] "S3"
Code
sloop::otype(unclass(lm0))
[1] "base"
Question

sloop exports functions s3_class() and s3_get_method()

  • Apply s3_class() to all members of lm0
  • What is the otype of autoplot() applied to an object of class lm?
Solution
Code
sloop::s3_class(lm0)
[1] "lm"
Code
sloop::s3_get_method(autoplot.lm)
function (object, which = c(1:3, 5), data = NULL, colour = "#444444", 
    size = NULL, linetype = NULL, alpha = NULL, fill = NULL, 
    shape = NULL, label = TRUE, label.label = ".label", label.colour = "#000000", 
    label.alpha = NULL, label.size = NULL, label.angle = NULL, 
    label.family = NULL, label.fontface = NULL, label.lineheight = NULL, 
    label.hjust = NULL, label.vjust = NULL, label.repel = FALSE, 
    label.n = 3, smooth.colour = "#0000FF", smooth.linetype = "solid", 
    ad.colour = "#888888", ad.linetype = "dashed", ad.size = 0.2, 
    nrow = NULL, ncol = NULL, ...) 
{
    p1 <- p2 <- p3 <- p4 <- p5 <- p6 <- NULL
    dropInf <- function(x, h) {
        if (any(isInf <- h >= 1)) {
            warning(gettextf("not plotting observations with leverage one:\n  %s", 
                paste(which(isInf), collapse = ", ")), call. = FALSE, 
                domain = NA)
            x[isInf] <- NaN
        }
        x
    }
    show <- rep(FALSE, 6)
    show[which] <- TRUE
    class(object$residuals) <- NULL
    if (is.null(data)) {
        plot.data <- ggplot2::fortify(object)
    }
    else {
        plot.data <- ggplot2::fortify(object, data = data)
    }
    n <- nrow(plot.data)
    plot.data$.index <- 1:n
    plot.data$.label <- rownames(plot.data)
    is_glm <- inherits(object, "glm")
    r <- residuals(object)
    w <- weights(object)
    if (any(show[2L:6L])) {
        s <- if (inherits(object, "rlm")) {
            object$s
        }
        else if (is_glm) {
            sqrt(summary(object)$dispersion)
        }
        else {
            sqrt(stats::deviance(object)/stats::df.residual(object))
        }
        hii <- stats::lm.influence(object, do.coef = FALSE)$hat
        r.hat <- range(hii, na.rm = TRUE)
        is_const_lev <- all(r.hat == 0) || all(diff(r.hat) < 
            1e-10 * mean(hii, na.rm = TRUE))
        fs <- dplyr::select_if(plot.data, function(x) is.character(x) | 
            is.factor(x))
        fs[[".label"]] <- NULL
        if (is_const_lev & ncol(fs) > 0) {
            plot.data$.nf <- stringr::str_wrap(interaction(fs, 
                sep = ":"), width = 10)
        }
        if (any(show[2L:3L])) {
            plot.data$.wresid <- if (is.null(w)) {
                r
            }
            else {
                sqrt(w) * r
            }
            plot.data$.wstdresid <- plot.data$.wresid/(s * sqrt(1 - 
                hii))
        }
        if (show[2L]) {
            ylim <- range(plot.data$.wstdresid, na.rm = TRUE)
            ylim[2L] <- ylim[2L] + diff(ylim) * 0.075
            qn <- stats::qqnorm(plot.data$.wstdresid, ylim = ylim, 
                plot.it = FALSE)
            plot.data$.qqx <- qn$x
            plot.data$.qqy <- qn$y
        }
    }
    label.fitted <- ifelse(is_glm, "Predicted values", "Fitted values")
    label.y23 <- ifelse(is_glm, "Std. deviance resid.", "Standardized residuals")
    if (is.logical(shape) && !shape) {
        if (missing(label)) {
            label <- TRUE
        }
        if (missing(label.n)) {
            label.n <- nrow(plot.data)
        }
    }
    plot.data <- flatten(plot.data)
    if (label.n > 0L) {
        if (show[1L]) {
            r.data <- dplyr::arrange(plot.data, dplyr::desc(abs(.resid)))
            r.data <- utils::head(r.data, label.n)
        }
        if (".wresid" %in% colnames(plot.data)) {
            wr.data <- dplyr::arrange(plot.data, dplyr::desc(abs(.wresid)))
            wr.data <- utils::head(wr.data, label.n)
        }
        if (any(show[4L:6L])) {
            cd.data <- dplyr::arrange(plot.data, dplyr::desc(abs(.cooksd)))
            cd.data <- utils::head(cd.data, label.n)
        }
    }
    .smooth <- function(x, y) {
        stats::lowess(x, y, f = 2/3, iter = 3)
    }
    .decorate.label <- function(p, data) {
        if (label & label.n > 0) {
            p <- plot_label(p = p, data = data, label = label, 
                label.label = label.label, label.colour = label.colour, 
                label.alpha = label.alpha, label.size = label.size, 
                label.angle = label.angle, label.family = label.family, 
                label.fontface = label.fontface, label.lineheight = label.lineheight, 
                label.hjust = label.hjust, label.vjust = label.vjust, 
                label.repel = label.repel)
        }
        p
    }
    .decorate.plot <- function(p, xlab = NULL, ylab = NULL, title = NULL) {
        p + ggplot2::xlab(xlab) + ggplot2::ylab(ylab) + ggplot2::ggtitle(title)
    }
    smoother_m <- ggplot2::aes_string(x = "x", y = "y")
    if (show[1L]) {
        t1 <- "Residuals vs Fitted"
        mapping <- ggplot2::aes_string(x = ".fitted", y = ".resid")
        smoother <- .smooth(plot.data$.fitted, plot.data$.resid)
        smoother <- as.data.frame(smoother)
        p1 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
        if (!is.logical(shape) || shape) {
            p1 <- p1 + geom_factory(geom_point, plot.data, colour = colour, 
                size = size, linetype = linetype, alpha = alpha, 
                fill = fill, shape = shape)
        }
        p1 <- p1 + ggplot2::geom_line(data = smoother, mapping = smoother_m, 
            colour = smooth.colour, linetype = smooth.linetype) + 
            ggplot2::geom_hline(yintercept = 0L, linetype = ad.linetype, 
                size = ad.size, colour = ad.colour)
        p1 <- .decorate.label(p1, r.data)
        p1 <- .decorate.plot(p1, xlab = label.fitted, ylab = "Residuals", 
            title = t1)
    }
    if (show[2L]) {
        t2 <- "Normal Q-Q"
        qprobs <- c(0.25, 0.75)
        qy <- stats::quantile(plot.data$.wstdresid, probs = qprobs, 
            names = FALSE, type = 7, na.rm = TRUE)
        qx <- stats::qnorm(qprobs)
        slope <- diff(qy)/diff(qx)
        int <- qy[1L] - slope * qx[1L]
        mapping <- ggplot2::aes_string(x = ".qqx", y = ".qqy")
        p2 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
        if (!is.logical(shape) || shape) {
            p2 <- p2 + geom_factory(geom_point, plot.data, colour = colour, 
                size = size, linetype = linetype, alpha = alpha, 
                fill = fill, shape = shape)
        }
        p2 <- p2 + ggplot2::geom_abline(intercept = int, slope = slope, 
            linetype = ad.linetype, size = ad.size, colour = ad.colour)
        p2 <- .decorate.label(p2, wr.data)
        p2 <- .decorate.plot(p2, xlab = "Theoretical Quantiles", 
            ylab = label.y23, title = t2)
    }
    if (show[3L]) {
        t3 <- "Scale-Location"
        mapping <- ggplot2::aes_string(x = ".fitted", y = "sqrt(abs(.wstdresid))")
        smoother <- .smooth(plot.data$.fitted, sqrt(abs(plot.data$.wstdresid)))
        smoother <- as.data.frame(smoother)
        p3 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
        if (!is.logical(shape) || shape) {
            p3 <- p3 + geom_factory(geom_point, plot.data, colour = colour, 
                size = size, linetype = linetype, alpha = alpha, 
                fill = fill, shape = shape)
        }
        p3 <- p3 + ggplot2::geom_line(data = smoother, mapping = smoother_m, 
            colour = smooth.colour, linetype = smooth.linetype)
        p3 <- .decorate.label(p3, wr.data)
        label.y3 <- ifelse(is_glm, expression(sqrt(abs(`Std. deviance resid.`))), 
            expression(sqrt(abs(`Standardized residuals`))))
        p3 <- .decorate.plot(p3, xlab = label.fitted, ylab = label.y3, 
            title = t3)
    }
    if (show[4L]) {
        t4 <- "Cook's distance"
        mapping <- ggplot2::aes_string(x = ".index", y = ".cooksd", 
            ymin = 0, ymax = ".cooksd")
        p4 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
        if (!is.logical(shape) || shape) {
            p4 <- p4 + geom_factory(geom_linerange, plot.data, 
                colour = colour, size = size, linetype = linetype, 
                alpha = alpha, fill = fill, shape = shape)
        }
        p4 <- .decorate.label(p4, cd.data)
        p4 <- .decorate.plot(p4, xlab = "Obs. Number", ylab = "Cook's distance", 
            title = t4)
    }
    if (show[5L]) {
        if (is_const_lev & ncol(fs) > 0) {
            t5 <- "Constant Leverage:\nResiduals vs Factor Levels"
            mapping <- ggplot2::aes_string(x = ".nf", y = ".stdresid")
            p5 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
            if (!is.logical(shape) || shape) {
                p5 <- p5 + geom_factory(geom_point, plot.data, 
                  colour = colour, size = size, linetype = linetype, 
                  alpha = alpha, fill = fill, shape = shape)
            }
            p5 <- p5 + ggplot2::geom_hline(yintercept = 0L, linetype = ad.linetype, 
                size = ad.size, colour = ad.colour) + ggplot2::expand_limits(x = 0)
            p5 <- .decorate.label(p5, cd.data)
            label.y5 <- ifelse(is_glm, "Std. Pearson resid.", 
                "Standardized Residuals")
            p5 <- .decorate.plot(p5, xlab = "Factor Level Combination", 
                ylab = label.y5, title = t5)
        }
        else {
            t5 <- "Residuals vs Leverage"
            mapping <- ggplot2::aes_string(x = ".hat", y = ".stdresid")
            smoother <- .smooth(plot.data$.hat, plot.data$.stdresid)
            smoother <- as.data.frame(smoother)
            p5 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
            if (!is.logical(shape) || shape) {
                p5 <- p5 + geom_factory(geom_point, plot.data, 
                  colour = colour, size = size, linetype = linetype, 
                  alpha = alpha, fill = fill, shape = shape)
            }
            p5 <- p5 + ggplot2::geom_line(data = smoother, mapping = smoother_m, 
                colour = smooth.colour, linetype = smooth.linetype) + 
                ggplot2::geom_hline(yintercept = 0L, linetype = ad.linetype, 
                  size = ad.size, colour = ad.colour) + ggplot2::expand_limits(x = 0)
            p5 <- .decorate.label(p5, cd.data)
            label.y5 <- ifelse(is_glm, "Std. Pearson resid.", 
                "Standardized Residuals")
            p5 <- .decorate.plot(p5, xlab = "Leverage", ylab = label.y5, 
                title = t5)
        }
    }
    if (show[6L]) {
        t6 <- "Cook's dist vs Leverage"
        mapping <- ggplot2::aes_string(x = ".hat", y = ".cooksd")
        smoother <- .smooth(plot.data$.hat, plot.data$.cooksd)
        smoother <- as.data.frame(smoother)
        p6 <- ggplot2::ggplot(data = plot.data, mapping = mapping)
        if (!is.logical(shape) || shape) {
            p6 <- p6 + geom_factory(geom_point, plot.data, colour = colour, 
                size = size, linetype = linetype, alpha = alpha, 
                fill = fill, shape = shape)
        }
        p6 <- p6 + ggplot2::geom_line(data = smoother, mapping = smoother_m, 
            colour = smooth.colour, linetype = smooth.linetype) + 
            ggplot2::expand_limits(x = 0, y = 0)
        p6 <- .decorate.label(p6, cd.data)
        p6 <- .decorate.plot(p6, xlab = "Leverage", ylab = "Cook's distance", 
            title = t6)
        g <- dropInf(hii/(1 - hii), hii)
        p <- length(stats::coef(object))
        bval <- pretty(sqrt(p * plot.data$.cooksd/g), 5)
        for (i in seq_along(bval)) {
            bi2 <- bval[i]^2
            p6 <- p6 + ggplot2::geom_abline(intercept = 0, slope = bi2, 
                linetype = ad.linetype, size = ad.size, colour = ad.colour)
        }
    }
    if (is.null(ncol)) {
        ncol <- 0
    }
    if (is.null(nrow)) {
        nrow <- 0
    }
    plot.list <- list(p1, p2, p3, p4, p5, p6)[which]
    new("ggmultiplot", plots = plot.list, nrow = nrow, ncol = ncol)
}
<bytecode: 0x609dfc00dcb8>
<environment: namespace:ggfortify>
Code
sloop::s3_get_method(augment.lm)
function (x, data = model.frame(x), newdata = NULL, se_fit = FALSE, 
    interval = c("none", "confidence", "prediction"), conf.level = 0.95, 
    ...) 
{
    warn_on_subclass(x, "augment")
    check_ellipses("level", "augment", "lm", ...)
    interval <- match.arg(interval)
    df <- augment_newdata(x, data, newdata, se_fit, interval, 
        level = conf.level)
    if (is.null(newdata)) {
        tryCatch({
            infl <- influence(x, do.coef = FALSE)
            df <- add_hat_sigma_cols(df, x, infl)
        }, error = data_error)
    }
    df
}
<bytecode: 0x609dfba87d00>
<environment: namespace:broom>
Question
Solution
Code
sloop::s3_methods_class("prcomp") |>
  gt::gt() |>
  gt::tab_header("Generics of class prcomp")
Generics of class prcomp
generic class visible source
augment prcomp FALSE registered S3method
autoplot prcomp FALSE registered S3method
biplot prcomp FALSE registered S3method
fortify prcomp FALSE registered S3method
plot prcomp FALSE registered S3method
predict prcomp FALSE registered S3method
print prcomp FALSE registered S3method
summary prcomp FALSE registered S3method
tidy prcomp FALSE registered S3method
Question
Solution
Code
sloop::s3_dispatch(autoplot(lm0))
=> autoplot.lm
 * autoplot.default
Question
Solution
Code
data(UCBAdmissions)
class(UCBAdmissions)
[1] "table"
Code
otype(UCBAdmissions)
[1] "S3"
Question
Solution
Code
sloop::s3_class(UCBAdmissions)
[1] "table"
Code
# sloop::s3_get_method(autoplot.table)
# sloop::s3_get_method(augment.table)

Programming with dplyr and ggplot2

We first aim at programming a function that takes as input a dataframe df, a column name col, and that, depending on the type of the column denoted by col, plots a histogram (for numerical column), a barplot (for factors), or raise an error of the column is neither categorical, nor numerical.

The function should return a ggplot object.

Here is a first attempt.

Let us first build a toy tibble.

Code
tb <- tibble( 
  col_num = rnorm(100), 
  col_fac = as_factor(sample(letters, 100, replace = T)), 
  col_ts = Sys.time() + duration(sample(1:20, 100, replace=T),units="days")
) 

tb |> 
  head()
# A tibble: 6 × 3
  col_num col_fac col_ts             
    <dbl> <fct>   <dttm>             
1  -0.425 z       2025-04-04 00:11:08
2   0.129 l       2025-03-23 23:11:08
3  -0.593 i       2025-03-23 23:11:08
4   0.554 s       2025-04-02 00:11:08
5  -1.21  e       2025-03-29 23:11:08
6   0.273 k       2025-03-25 23:11:08
Solution
Code
gg_obj <-  function(df, col){
  
  vct <- df[[col]]
  tp <- class(vct)

  if (tp != "numeric" & tp !="factor") {
    stop(paste0(col, " is of wrong type!"))
  }

  p <- ggplot(df) + 
    aes(x=.data[[col]])

  if (tp=="numeric") {
    p <- p + geom_histogram()
  } else {
    p <- p + geom_bar()
  }

  p  
}
1
List component is accessed by name. col is a string.
2
.data is a pronoun for the dataframe component of the ggplot object
Code
(
  gg_obj(tb, "col_num") +
    labs(
      title= "Histogram", 
      subtitle= "Numerical column")
) + (
  gg_obj(tb, "col_fac") +
    labs(
      title= "Barplot",
      subtitle= "Factor column"
    )
)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Question
  • Pass more optional arguments to geom_... (use ellipsis ...)
  • Avoid quoting the column name
Code
gg_obj_2 <-  function(df, col, ...){
  params <- list(...)
  vct <- pull(df, {{col}})
  tp <- class(vct)[1]

  if (tp != "numeric" & tp !="factor") {
    stop("column is of wrong type!")
    return
  }

  p <- ggplot(df) + 
    aes(x={{col}})

  if (tp=="numeric") {
    p <- p + geom_histogram(...)
  } else {
    p <- p + geom_bar(...)
  }

  p  
}
1
standard technique for programming with tidy-evaluation: embracing.
2
using the ellipsis
Code
(
  gg_obj(tb, "col_num") +
    labs(
      title= "Histogram", 
      subtitle= "Numerical column")
) +
(
  gg_obj_2(tb, col_num, fill="red") +
    labs(
      title= "Histogram, tuning color", 
      subtitle= "Numerical column")
)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Question
Question

How could you add a geom_point() layer to each element of the following list?

Code
plots <- list(
  ggplot(mpg, aes(displ, hwy)),
  ggplot(diamonds, aes(carat, price)),
  ggplot(faithfuld, aes(waiting, eruptions, size = density))
)

From R Advanced Programming

Solution
Code
plots |> 
  map(\(x) x + geom_point()) |>
  patchwork::wrap_plots()

Inside lm()

Question

In classes like lm, prcomp, … we have a member called call. What does it represent? How is it constructed?

First, read the code of lm.

> lm 
function (formula, data, subset, weights, na.action, method = "qr", 
    model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, 
    contrasts = NULL, offset, ...) 
{
    ret.x <- x
    ret.y <- y
    cl <- match.call()
    mf <- match.call(expand.dots = FALSE)
    m <- match(c("formula", "data", "subset", "weights", "na.action", 
        "offset"), names(mf), 0L)
    mf <- mf[c(1L, m)]
    mf$drop.unused.levels <- TRUE
    mf[[1L]] <- quote(stats::model.frame)
    mf <- eval(mf, parent.frame())
    if (method == "model.frame") 
        return(mf)
    else if (method != "qr") 
        warning(gettextf("method = '%s' is not supported. Using 'qr'", 
            method), domain = NA)
    mt <- attr(mf, "terms")
    y <- model.response(mf, "numeric")
    w <- as.vector(model.weights(mf))
    if (!is.null(w) && !is.numeric(w)) 
        stop("'weights' must be a numeric vector")
    offset <- model.offset(mf)
    mlm <- is.matrix(y)
    ny <- if (mlm) 
        nrow(y)
    else length(y)
    if (!is.null(offset)) {
        if (!mlm) 
            offset <- as.vector(offset)
        if (NROW(offset) != ny) 
            stop(gettextf("number of offsets is %d, should equal %d (number of observations)", 
                NROW(offset), ny), domain = NA)
    }
    if (is.empty.model(mt)) {
        x <- NULL
        z <- list(coefficients = if (mlm) matrix(NA_real_, 0, 
            ncol(y)) else numeric(), 
                  residuals = y, 
                  fitted.values = 0 * y, 
                  weights = w, 
                  rank = 0L, 
                  df.residual = if (!is.null(w)) sum(w != 0) else ny
              )
        if (!is.null(offset)) {
            z$fitted.values <- offset
            z$residuals <- y - offset
        }
    }
    else {
        x <- model.matrix(mt, mf, contrasts)
        z <- if (is.null(w)) 
            lm.fit(x, y, offset = offset, singular.ok = singular.ok, 
                ...)
        else lm.wfit(x, y, w, offset = offset, singular.ok = singular.ok, 
            ...)
    }
    class(z) <- c(if (mlm) "mlm", "lm")
    z$na.action <- attr(mf, "na.action")
    z$offset <- offset
    z$contrasts <- attr(x, "contrasts")
    z$xlevels <- .getXlevels(mt, mf)
    z$call <- cl
    z$terms <- mt
    if (model) 
        z$model <- mf
    if (ret.x) 
        z$x <- x
    if (ret.y) 
        z$y <- y
    if (!qr) 
        z$qr <- NULL
    z
}
<bytecode: 0x55564224e930>
<environment: namespace:stats>
Question

Have a look at function match.call()

Code
e <- match.call(get, call("get", "abc", i = FALSE, p = 3))

is_expression(e)
[1] TRUE
Code
is_call(e)
[1] TRUE
Code
lobstr::ast(match.call(get, call("get", "abc", i = FALSE, p = 3)))
█─match.call 
├─get 
└─█─call 
  ├─"get" 
  ├─"abc" 
  ├─i = FALSE 
  └─p = 3 

Let us explore the next toy example.

Code
fun <- function(x, lower = 0, upper = 1) {
  structure((x - lower) / (upper - lower), CALL = match.call())
}


w <- fun(4 * atan(1), u = pi)

type_of(w)
Warning: `type_of()` is deprecated as of rlang 0.4.0.
Please use `typeof()` or your own version instead.
This warning is displayed once every 8 hours.
[1] "double"
Code
w
[1] 1
attr(,"CALL")
fun(x = 4 * atan(1), upper = pi)
Code
str(w)
 num 1
 - attr(*, "CALL")= language fun(x = 4 * atan(1), upper = pi)
Code
x <- attr(w, 'CALL')

typeof(x)
[1] "language"
Code
is_call(x)
[1] TRUE
Code
eval(x)  
[1] 1
attr(,"CALL")
fun(x = 4 * atan(1), upper = pi)
Code
# try also eval_tidy(), eval_bare

Data masking and environments

Question
Solution

Tidy evaluation

Question

What is quasi-quotation?

Keep the rlang cheatsheet around.

Question

Explain the difference between an expression and a quosure

Question

References

Programming with ggplot